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ABSTRACT 


Analysis of atmospheric motion over obstacles on plane surfaces 
is carried out in the present study to compute simulated wind fields 
over terrain features. Emphasis is on a semielliptical, two-dimensional 
geometry. Numerical simulation of flow over rectangular geometries is 
also discussed. The numerical procedure utilizes a two-equation turbu- 
lence model^ and the selection of the necessary constant coefficients in 
the model is considered an important part of the present study. 

In the present approach the partial differential equations for 
the vorticity, stream function, turbulence kinetic energy, and turbu- 
lence length scale are solved by a finite-difference technique. The 
numerical solutions have been compared with available experimental data; 
agreement is good. 

It is found that the mechanism of flow separation induced by a 
semi ell ipse is the same as in the case of flow over a gradually sloping 
surface for which the flow separation is caused by the interaction 
between the viscous force, the pressure force, and the turbulence level. 
For flow over bluff bodies, e.g., solid fences, a downstream 
recirculation bubble is created due to the inability of the flow to 
negotiate with the abrupt change of the surface shape. Increasing the 
aspect ratio and/or increasing the turbulence level results in flow 
reattachment close behind the obstacle. 
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1.0 INTRODUCTION 


Analysis of atmospheric motion over obstacles on plane surfaces 
is carried out in the present study to compute simulated wind fields 
over terrain features. Emphasis is on a semielliptical, two-dimensional 
geometry. Numerical simulation of flow over rectangular geometries is 
also discussed. The numerical procedure utilizes a two-equation turbu- 
lence model and the selection of the necessary constant coefficients in 
the model is considered an important part of the present study. 

Although the terrain features analyzed are highly idealized 
geometries which are not exactly duplicated by natural terrain, the 
results present are sufficiently representative of the real world to 
provide insight into the various hazards which are associated with 
aircraft flight in the vicinity of both natural and man-made objects. 
Controlled flight over terrain, such as landing off the runway environ- 
ment or flying over rising terrain, continues to be a safety issue. 

In recent accidents where terrain was a major factor, the air traffic 
control system on the aircraft was equipped with the capabilities 
that should have prevented the accident. This suggests that there 
were certain unknown factors which may include complex and unexpected 
wind patterns which for this safety issue need re-examination. 

Moreover, operation of low-speed helicopters and V/STOL aircraft in 
large metropolitan areas is also affected by zones of recirculating 
wind fields, regions of large velocity fluctuations, and fields of 
vorticity induced by the atmospheric motion near buildings. 
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These complex wind fields can be extremely hazardous during landing and 
takeoff. 

There are many other technological developments for which 
analytical solutions of flow over terrain features have practical appli- 
cations. The dispersion of air pollution is significantly influenced by 
atmospheric motion. In complex terrain, the movement and deposition of 
the pollutant can be significantly altered by a hill or building in the 
vicinity of the stack release. The siting of wind turbine generators 
requires knowledge of regions of high wind speed and avoidance of 
regions where recirculation or highly turbulent flow may occur. Addi- 
tionally, the wind load on structures and wind erosion can be increased 
or decreased by the wind characteristics in the vicinity of surface 
irregularities. These fields of technology also require an understanding 
of the atmospheric flow properties about surface obstacles. 

It is well known that surface wind motion is markedly affected 
by the abrupt change of terrain features [1].^ The phenomena of flow 
acceleration, of flow separation, and of augmented turbulence in the 
wind field are observed effects. The effect of shear in the approaching 
flow and the strong pressure gradients near the barriers create zones 
of recirculation both upstream and downstream in the vicinity of the 
obstacle. High turbulence is generated in the shear layer bordering the 
recirculation regions due to the interaction between the turbulence 
stresses and the deformation of the mean flow field. Momentum in the 
shear layer diffuses into the wake and into the quasi-potential flow 

^Numbers in brackets refer to similarly numbered references in 
the Bibliography. 
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outside the wake, setting the wake fluid into motion and smoothing out 
the sharp velocity discontinuities. Downstream of the obstacles, the 
diffusion of momentum gradually thickens the shear layer until the inner 
flow is blended with the outer flow, forming a new boundary layer. 

The nature of these flow fields has been studied experimentally 
as well as analytically. Reported results of some of these studies are 
reviewed in Section 2.0. Due to the complexity of the flow characteris- 
tics near the obstacles, it is necessary to solve the complete equations 
of motion if a realistic recirculation flow field is to be computed. 
Practical methods, however, for solving such equations require the use 
of numerical modeling techniques. 

The numerical procedure employed in this paper is an elliptical 
finite-difference algorithm initially developed by Gosman, et al. [2]. 

The Navier-Stokes equations have been expressed in terms of stream func- 
tion and vorticity. In order to close the governing equations, the 
Reynolds shear stresses appearing in the time-averaged equations of 
motion are expressed by the product of a scalar eddy viscosity and a 
strain rate. The eddy viscosity is expressed as a function of the turbu- 
lence kinetic energy and the turbulence length scale. These quantities 
are determined, respectively, from two corresponding transport equations 
(two-equation turbulence model [3]). The computation of turbulent flows 
based on current turbulence models in all cases requires the determina- 
tion of one or more constant coefficients. One goal of this study is to 
derive the relationships of these constants from the knowledge of turbu- 
lence characteristics available at the present time. Details of the 
turbulence model employed in this study are discussed in Section 3.0. 
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Solving partial differential equations by numerical techniques 
requires accurate numerical presentation of boundary conditions. When 
the boundary is coincident with a coordinate line, the finite-difference 
expression at and adjacent to the boundary can be applied at grid 
points constructed on the coordinates without the need of interpolation 
between grid points. In this study, the governing differential equa- 
tions are written in an orthogonal curvilinear coordinate system. The 
curved coordinate surface can then be fit to most terrain shapes of 
interest. Although this increases the size of the governing equations 
slightly, little additional complication is introduced into the calcula- 
tion procedure. The governing equations as well as the coordinate system 
used in this study are described in Section 4.0. The numerical proce- 
dures and the imposed boundary conditions are given in Section 5.0. 

The computation has been carried out for atmospheric flow over 
two-dimensional semiellipses, due to the resemblance of their geometries 
to the shapes of hills. An analysis of the effect of the geometry of 
the obstacles on the wind field is carried out in the present study by 
varying the aspect ratio of the semiellipse. The atmospheric flow is 
considered neutrally stable. The concept of a neutrally stable atmo- 
sphere is meaningful in strong wind conditions [1,4], which is of primary 
importance in this study. Good agreement between the numerical predic- 
tions and the experimental data obtained from both atmospheric boundary 
layer wind tunnels [5,6,7] and full-scale field measurements [8,9] is 
shown in Section 6.0. The comparisons are carried out for flows over 
solid fences, rectangular blocks and escarpments, in addition to the 
case of flow over semielliptical bodies. 
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Parameters influencing the flow field in the vicinity of surface 
obstacles are also discussed in Section 6.0. It is found that the 
primary flow parameters affecting the flow field are surface roughness 
scale, the shape of the obstacles as well as the undisturbed upstream 
turbulence level. Increasing the scale length of the surface roughness 
results in reduced wind speed in the lower layer of the atmosphere due 
to the increase in surface frictional drag. Thus, the streamline is dis- 
placed to a greater height when the flow passes over the obstacle with 
larger surface roughness. For the flow conditions investigated in the 
present study, flow separation is not obtained for flow over semi ellipses 

5 

when the molecular Reynolds number, Re, is larger than 10 . Conditions 
with Re = 100, however, show that separated flow is induced by semi- 
elliptical obstacles. The mechanism of flow separation induced by a 
semiellipse is the same as in the case of flow over a gradually slop- 
ing surface for which the flow separation is caused by the interaction 
between the viscous force, the pressure force, and the turbulence level. 
For flow over bluff bodies, e.g., solid fences and rectangular blocks, 
a large downstream recirculation bubble is created due to the inability 
of the flow to negotiate with the abrupt change of the surface shape. 

It is found that the aspect ratio of the block and the turbulence level 
of the approaching flow cause significant effect on the size of the 
downstream recirculating flow region. Increasing the aspect ratio 
and/or increasing the turbulence level result in flow reattachment 
close behind the obstacle. Details of flow separation for semi- 
ellipses as well as for bluff bodies are discussed in Section 6.5. 
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2.0 REVIEW OF PREVIOUS WORKS 


The purpose of this section is to review the existing theoret- 
ical analyses and experimental results concerning the aerodynamics of 
surface obstructions in turbulent flow fields. Frost [10], has compiled 
an extensive survey of flow properties around man-made obstacles to the 
wind. The obstacle exerts a drag force on the wind field and distorts 
the flow approaching and passing over it. Recirculation regions are 
created due to the strong pressure gradients near the obstacle. The 
separation and reattachment points in two-dimensional laminar flow fields 
correspond to the location where the normal velocity gradient is zero; 
thus the wall shear is zero. In a turbulent flow field, this is also 
assumed to be true in terms of the ensemble time-averaged quantities [11]. 

Various approaches have been used to theoretically describe 
certain features of turbulent flow over barriers. Tani [12] used a 
diffusion equation to investigate the velocity distributions behind a 
permeable fence. Velocities in the vicinity of the fence were not 
accurately predicted. Plate [13] and Chang [14] indicate that the 
velocity distributions in the shear layer bordering the recirculation 
bubble behind an obstacle agree with prediction from mixing-layer 
theory. No flow separation, however, is taken into account. The 
agreement may be due to the arbitrary adjustment of coefficients used 
in the theory. 

Counihan, et al . [15] apply a small perturbation technique to 
study flow over a low block. Their conclusion that similarity prevails 
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and that the velocity deficit behind the block correlates with a single 
curve is also supported by their experimental work. Some evidence [5], 
however, is given that such similarity in general may not exist. It is 
unlikely that details of the flow properties in the recirculation 
regions near the barrier can be accurately predicted by using small 
perturbation techniques. 

Several investigators, e.g., Taylor [16], Jackson and Hunt [17], 
and Deaves [18], are studying turbulent flow over gently rolling, low 
hills. The application of their solution to large or steep hills must 
be treated with caution. Frost, et al . [19] extended a turbulent boundary- 
layer theory to approximate the atmospheric motion over semiel 1 iptical 
hills. The Reynolds stresses are assumed to be important only close to 
the ground, and the pressure distributions driving the boundary layer 
are estimated from inviscid flow theory. This assumption is supported 
by the work of Jackson and Hunt [17] for flow over low hills. 

There has been a fair amount of work recently on the problem 
of flow over circular or elliptical cylinders, most of which assumes low 
Reynolds number flow, e.g., Lin and Lee [20] and Haussling [21], or low 
turbulence high Reynolds number flow, e.g., Guven, et al . [22]. None of 
these studies are easily extended to the analysis of turbulent flow 
over barriers. 

At present, most understanding of turbulent flow about surface 
obstructions is obtained from empirical observations. Frost and Shieh 
[1] have made an extensive survey of the available literature. The 
following summarizes their survey. 
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The presence of a bluff body in the flow field retards the 
approaching wind and creates a separation bubble in front of the 
obstacle. The pressure on the front face of the barrier is strong due 
to the conservation of energy. Outside the upstream separation bubble 
these strong pressure gradients displace the flow and accelerate it 
as it passes over the obstacle. A lee separation bubble is created 
if the flow cannot negotiate the change in slope of the barrier. Mass 
flow into and out of the recirculation zone is conserved. At the 
reattachment point the flow splits, part flowing into the bubble and 
part flowing downstream. An equal amount of mass is in turn entrained 
into the free stream from the bubble through the turbulent transport 
processes in the shear layer bordering the recirculation zone [23]. 
Downstream of the reattachment point, a new boundary layer develops. 

The retarded wind profile readjusts to the local surface conditions, 
thus recovering to a new, fully developed profile at a distance far 
downstream. The behavior of this profile has been studied with 
boundary-layer theory [24]. 

The flow properties near an obstacle depend on the geometry of 
the obstacle, the local surface conditions, and the nature of the 
upstream wind. Behind the obstacle, a closed recirculating flow region 
(recirculation bubble) bounded by a stagnation streamline is generally 
accepted as a description of the separation zone in two-dimensional 
flow. The concept of a closed recirculation bubble can also be 
applied to the case of thick boundary layers, e.g., atmospheric 
boundary layer, if the turbulent transport process is more dominant 
than the mean convective transport process [25]. Strong turbulence 
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in the undisturbed flow causes a smaller recirculation zone on the lee 
side of most obstacles [26]. 

One of the characteristics of the wake generated by an obstacle 
is the velocity deficit in the wake. The decay rate of the mean 
velocity deficit is independent of the obstacle size, but the magnitude 
of the disturbance is greater for larger obstacles; hence the wake 
persists farther downstream [11]. 

The shape of a barrier has a definite effect on the lee side 
velocity distributions. A solid vertical sharp-edged fence generally 
creates a stronger disturbance than a porous fence, although the dis- 
turbance may be carried farther downstream behind a porous fence. 
Acceleration of the flow occurs in a region located a distance two to 
three fence heights downstream and slightly above the lee side recir- 
culation bubble. Also, large reductions in velocity occur close to the 
ground immediately behind the fence for solid vertical plate fences. 
Increasing the permeability of the fence causes less reduction of the 
flow behind the fence as well as weaker wind speed in the region of 
flow acceleration over the bubble. The area of reduced wind speed 
behind a porous fence persists greater distances downstream due to more 
momentum being convected into the leeward area through the fence 
openings. The windward flow patterns, however, are insignificantly 
affected by the barrier porosity. 

The maximum wind acceleration expressed in terms of the ratio 
of the local wind speed relative to the undisturbed value measured at 
the same height above the local surface occurs at the crest of a 
surface feature in the absence of flow separation [7]. When flow 
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separation occurs, the location of the maximum speed lies a small dis- 
tance downstream of the barrier and slightly above the separation region, 
depending on the barrier shape, the local surface conditions and the 
turbulence of the approaching flow. 

Vlhen the wind approaches a three-dimensional feature, flow 
acceleration occurs not only on top of the feature but also in regions 
near the sides. In a stably stratified atmosphere, the wind tends to 
pass by an obstacle rather than passing over it. Due to this effect, 
the velocity on top of the obstacle is significantly reduced for flow 
over terrain features with small aspect ratios. However, for flow over 
a long ridge the velocity at the center portion of the ridge is not 
significantly affected by atmospheric stratifications due to quasi - 
two-dimensional flow conditions generally prevalent at that portion 
[7]. Additionally, the flow separation point is observed to move down- 
wind under stable conditions [1]. 

High turbulence generally occurs in a shear layer originating 
from the sharp edge of a bluff body. This high turbulence diffuses 
vertically as flow is convected downstream, thus thickening the shear 
layer. The velocity fluctuation in the streamwise direction is not as 
sensitive to the barrier shapes as is the mean velocity [27]. High 
turbulence intensities, because of their definition, occur in regions 
of low wind speed, normally lying close to the rear face of an obstacle. 
Absolute turbulence, however, is generally higher in the regions of high 
wind speed. In the center of a long fence, where quasi-two-dimensional 
conditions are achieved, the streamwise fluctuation decreases with 
the increasing angle between the wind direction and the fence. 


10 



In general, the magnitudes of the vertical fluctuations are as large 
as the lateral fluctuations [27]. 
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3.0 TURBULENCE MODELS 


3.1 Mean Flow Equations 

In the study of atmospheric flow over surface obstacles, the 
flow associated with separation and reattachment is so complicated that 
the general equations of fluid motion, the Navier-Stokes equations, 
must be solved if realistic results are to be obtained. Atmospheric 
flows are generally treated as incompressible turbulent flow, i.e., the 
fluctuation of fluid density is negligible in comparison with the mean 
density in the statistical steady state. This is especially true when 
the atmosphere is neutrally stable. Based on the assumption of 
incompressible flow, the mean velocity and pressure are governed by 
the equations: 

Continuity equation: 



= 0 


(3.1) 


Momentum equation: 


8V. 9V. 

— I + V . — - 
9t J 9x- 
\) 


1 ^ 

p 9x^. 9Xj 


V 



V .V . 
1 J 


(3.2) 


where q. is the external force acting on the fluid element in the i 
direction and v is the kinematic viscosity of the fluid. 
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For a stable or an unstable atmosphere where thermal effects are sig- 
nificant, the influence of density variations due to temperature changes 
is often approximated with the Boussinesq approximation. The density 
effect thus appears only in the force term of Equation 3.2. The 
temperature distribution is then governed by a scalar transport equation 
for turbulent flow [28]: 

Scalar equation: 


IE 

at 


+ 



i 

p 3Xj 


(-pv^f) 


+ s, 


(3.3) 


where F is a scalar function, e.g., temperature; f is the fluctuation 
component of F; and Sp is a volumetric source term. Since the effect 
of heat transfer is negligible if the wind is strong [4], which is 
assumed in this study, the thermal effects on the flow are neglected, 
and therefore, the governing equations in this study are Equation 3d 
3o2. The external force term of Equation 3.2 is neglected due to its 
small effect compared with the shear stresses. 


3.2 Closure Problem 

The main problem in solving Equations 3.1 and 3.2 and thus cal- 
culating turbulent flows is the determination of the Reynolds stresses, 
-pv.v-. One can derive an exact transport equation for the Reynolds 

* w 

stresses from the Navier-Stokes equations. This derivation, however, 

simply introduces a third-order correlation term, v^.v^Vj^, [28], which 

still does not permit the system of equations to be closed. Continuing 
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the process by deriving exact equations for higher order velocity 
correlations results in still higher order correlations, and thus the 
set of equations can be closed. Thus, at some point a turbulence model 
must be introduced which approximates the turbulent correlations at 
some order in the hierarchy as functions of the lower order correlations 
and/or mean flow quantities. Turbulence is thus simulated by a turbu- 
lence model ranging from a simple algebraic model of mixing length to a 
more complicated differential transport equation model for the higher 
order correlations [3]. The additional equations, differential and/or 
algebraic equations, used for turbulence modeling together with Equations 
3.1 and 3.2 form a closed set of equations. The turbulence model used 
in this study is described in the following. 


3.3 Concept of Eddy Viscosity 

The eddy viscosity concept of Boussinesq [23] is adapted to 
express the Reynolds stresses as: 


-pv.v. 


pv, 


3V. 


3V.) 

-|- sJ- 

3X. 

j y 



(3.4) 


where 6.. is the Kronecker delta (6.j = 1 if i = j, otherwise S.. = 0), 

1 J I J I J 

and K is the kinetic energy of turbulence, l/2(v^. v . ) . The eddy vis- 
cosity, v^, used in Equation 3.4 is based on the assumption of an 
analogy between momentum transport by turbulence and momentum transport 
by molecular motion. Although this concept has frequently been criti- 
cized as physically unsound, it is found to work well in several flows 
[29]. 
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The concept of eddy viscosity has been used in isotropic turbu- 
lent flows by assuming a scalar quantity for v^. For nonisotropic 
turbulence fields, different eddy viscosities are sometimes used for 
turbulent transport of momentum in different directions [29]. In certain 
flows, e.g., wall jets and asymmetric wall shear layers, negative 
values of occur where the Reynolds stresses have the opposite sign 
to that of the velocity gradients (see Equation 3.4). This, of course, 
is not physically meaningful. In most cases, however, the regions with 
negative eddy viscosity are small, and therefore the error induced in 
those areas is of little practical importance [3]. The application of 
Equation 3.4, however, simply shifts the requirement of Reynolds 
stress modeling to the problem of eddy viscosity modeling. 

To overcome this problem, is expressed as being proportional 
to a velocity scale, /K, and a length scale, L, characterizing the 
turbulent motions. Prandtl and Kolmogorov [3] propose: 

V. = C /K L (3.5) 

t y 

where is a constant. The kinetic energy of turbulence, K, and its 
length scale, L, are computed from transport equations. This model, 
called a two-equation model, is employed in this study and is discussed 
in the next section. 

3 . 4 T he Two-Equation Mod els 

Launder and Spalding [3] indicate that a variable A constructed 
by any combination of K and L, i.e., A = where m and n are real 
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numbers, can be used in the transport equations for the two-equation 
models. Rodi and Karlsruhe [29], however, suggests that a transport 
equation for L should not be used due to the difficulties associated 
with deriving a generally valid formula. References [3,30,31] counter 
this and indicate that a common form for the L transport equation does 
exist. A transport equation for L having a general form like Equation 
3.3 has been used in References [2,11,32,33]. Evidence that realistic 
results can be achieved by using a K-L model for complicated separating 
and reattaching flow phenomena is reported in [11,32]. One subject of 
this study is to demonstrate the validity of the K-L model for calcul- 
ating turbulent recirculating flows. 

Forms of the two-equation model other than that used here are 
reported [3]. A common characteristic of all two-equation models is 
that a set of constant coefficients which appears in the transport equa- 
tions must be determined empirically or by inductive reasoning. One 
semi -empirical procedure for determining the constants assumed the 
turbulence is always in local equilibrium [3,32]. This approach, how- 
ever, can break down in regions where the turbulence production rate 
dominates the dissipation rate. This study proposes an alternate 
method to determine these constants based on presently available 
knowledge of turbulent motions. This method is discussed in detail 
in Section 3.6. 

3.5 K-L Two-Equation Model 

Derivations for the transport equation for K and for L are given 
in Hinze [28] and Mellor and Herring [30]. A brief description of their 
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development is given here. In the flow field, the convective and diffu- 
sive transport of turbulence, i.e., the history of the turbulence, is 
accounted for by solving two differential transport equations. These 
are the transport equations for turbulence velocity scale and for 
turbulence length scale. The physically most meaningful velocity scale 
is /K, where the turbulence kinetic energy is defined as half of the 
sum of the mean square values of the velocity fluctuation components, 
(l/2)v^v^. The transport equation for K can be derived from the Navier- 
Stokes equations [28]: 
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(3.6) 


I - = Diffusion 

II = Pi, = Production 

l\ 

III = Viscous transport 
IV = e,, = Dissipation 


The viscous transport term in Fquation 3.6 can be written as 
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since the terms 


* si 


and 


9^v . 

V . ^ 

1 3x.9x, 


1 


are zero for incompressible fluid. Since the viscous term is generally 
small except when spatial gradients of the Reynolds stress components 
are extremely large, it is generally neglected. The diffusion term in 
Equation 3.6 resulting from the interacting of velocity and pressure 
fluctuations is assumed to be subject to gradient diffusion. The 
diffusion term is then written 


D 


K 



^ IL 

0 ^. 3x. 

I ^ U 


(3.7) 


where the Schmidt number is a constant which is discussed in 
Section 3.6. Since the energy cascade process is independent of 
the molecular viscosity except at the final stage where the turbulence 
kinetic energy is converted into heat by viscous dissipation, the 
dissipated turbulence energy is modeled by 




(3.8) 


where is a constant which is discussed in Section 3.6 . Equation 
3.8 involves Kolmogorov's assumption [31] that the dissipation rate 
of K is determined by the energy-containing motion at high Reynolds 
number. The length scale, L, characterizes the size of the large, 
energy-containing eddies of the turbulence field. 
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The K equation based on the above assumptions is thus modeled as 
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(3.9) 


The length scale characterizing the size of energy-containing 
eddies, in turn, is assumed to be subject to transport processes and is 
also determined from a differential transport equation. As mentioned 
in the previous section, the L equation has the general form of equation 
3.3. The rate of change of L is balanced by the convective transport, 
the diffusive transport, D^, and the production and dissipation rate of 
L in the energy cascade process, i.e., and respectively, hence 


+ V = n + p + p 
9t j 9x. 


(3.10) 


It is assumed that the diffusive transport of L is also gradient-driven. 
The D|^ term is therefore modeled as 



t 9L 




(3.11) 


where a[_ is a constant which will be discussed in Section 3.6. 
Considering that increasing the turbulence kinetic energy dissipation 
increases the fatality rate for small eddies and thus effectively 
increases the eddy size, the production term for L, Pj^, should be 
related to the dissipation of turbulence kinetic and some characteristic 


time scale. Hence, it is argued that 
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which reduces to 



(3.12) 


where is a constant which will be discussed in Section 3.6. 

The vortex stretching connected with the energy cascade reduces the 
eddy size; thus, the rate of dissipation of L, is subject 
to the interaction of the mean motion and the turbulence fluctua- 
tions. 




t 3x. 

sJ . 


L K 


-1 


(3.13) 


where Cp is a constant which will be discussed in Section 3.6. 
With Equation 3.5 and Equations 3.11 through 3.13, Equation 3.10 
becomes 


~ + V 
9t j 9x. 
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(3.14) 


The length scale equation. Equation 3.14, traces the historical develop- 
ment of the length scale at a given location and therefore is strongly 
dependent on the initial value of L. 

Both Equations 3.9 and 3.14 involve several constants which must 
be determined before these equations can be solved. The determination 

of these constants is discussed in the next section, 
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3 . 6 Constants of the K-L Model 


The K-L model described in the previous section contains six 
coefficients. They are o^, C^, and Conceptually, these 

coefficients are functions of dimensionless parameters, e.g., Reynolds 
number, which govern the motion of the flow [2]. In high Reynolds 
number flows, however, the coefficients are generally treated as con- 
stants. For fully developed turbulent flows, the Schmidt numbers 
and 0 ^ used in Equations 3.9 and 3.14, respectively, are generally 
assumed to be of the order of unity [2]. 

a,_ = 1.0 (3.15) 

Indirect support of this assumption is given by Reynolds [34], who 
shows that the ratio of the eddy viscosity to the eddy diffusi vi ties 
of enthalpy and of other properties such as fluid contaminants is close 
to unity for shear flows. 

In the analysis of turbulent flows using two-equations models, 
determination of the constants C^, Cq, C^ and Cj^ depends either on 
experiment and/or on an exact analysis of some simple flow. Launder and 
Spalding [3] proposed a wall region approach in which the turbulence 
convection and diffusion are negligible. They analytically determine 
the constant coefficients used in their two-equation model in this 
region. Frost, et al . [32] also used this approach in determining the 
coefficient used in their K-L model. Moreover, the Prandtl mixing 
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length, L = kZq, was specified as a boundary condition for the transport 
equation of length scale. Equation 3.14. Realistic results are obtained 
by applying this model to calculate atmospheric flow over forward-facing 
: steps [32]. For flow over rectangular blocks, the wall region approach 
gives unrealistic results in both mean wind field and turbulence field 
[11]. Shi eh, et al. [11] investigated the influence of the coefficients, 
Cy, field computed and then proposed 

= 0.416, Cp = 0.416, = 1.44 and = 0.6 for calculating flows 

over a rectangular block. Realistic results in both mean wind field 
and turbulence field were obtained [11 ,26]. However, validation of 
the coefficients proposed by Shi eh, et al . [11] in other flow situa- 
tions requires further investigation. 

In the present study, the constant coefficients C , Cn» Cr and 

put 

C|^ described in the following are based on an understanding of the 
properties of turbulence as influenced by their transport process. 

Thus the coefficients are expected to have general applicability. The 
coefficients developed in this study are successfully applied to the 
cases of flow over semielliptical obstacles, fences, rectangular blocks 
and escarpments. Good agreement between the numerical solutions 
and the empirical data of both wind tunnel and full-scale field measure- 
ments are obtained. Comparisons of experiments and analyses are shown 
in Section 6.0. The coefficients are developed as follows. 

Taking the ratio of the production terms to the dissipation 
term in Equations 3.9 and 3.14, one obtains 
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( 3 . 17 ) 
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To establish some physical explanation of the relationship between the 
turbulence kinetic energy and the turbulence length scale, one can 
argue that the growth of the turbulence length scale is a result of the 
dissipation of energy, which creates the biggest fatality rate for 
small eddies. Thus, when dissipation exceeds production, the length 
scale, which is taken to represent some effective mean of the collection 
of eddy sizes, will increase with the depletion of the smaller eddies. 
The reduction in length scale, in turn, may be thought of as being 
caused by the tendency of shear stress to rupture the larger eddies. 
Therefore, when the production of turbulence kinetic energy exceeds 
dissipation, the length scale is expected to decrease, and when dissipa- 
tion of turbulence kinetic energy exceeds production the length scale 
is expected to grow. With this in mind it is proposed that 



which gives the relationship 



1.0 


( 3 . 18 ) 


This assumption is obviously valid for the local equilibrium case, i.e., 
Pk = £|^ and P|^ e However, it is not restricted to only equilibrium 
flows. 
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Due to the interaction between the mean motion and the turbulent 
motion, energy is extracted from the mean flow and converted into heat. 
The local production of turbulence energy, however, is not always 
balanced by the local viscous dissipation; thus the turbulence does not 
necessarily gain energy from the mean flow at the same rate as it loses 
energy through viscous dissipation. 

Conceptually, the ratio (C^/Cq) appearing in Equation 3.16 can 
be determined from the ratio and the local conditions of the 

turbulence and of the mean flow. Rodi [29] proposed an equation for 
evaluating the effect of the ratio on the coefficients used in 

a K-G|^ model. His coefficients vary significantly with the distortion 
rate of turbulence. However, an asymptotic constant value of is 
achieved if P|^/C]^ > 1*0. The distributions of K and L are shown [11] 
to be strongly dependent on the ratios C^/Cq and respectively, 

rather than on the magnitude of the constants themselves. However, 
close behind a rectangular block, x = 1.0 H, where the mean flow is 
highly distorted, the distributions of K and L are rather insensitive 
to the value of these ratios (Figure 3.1a). On the other hand, in the 
relaxation flow region far downstream from the block, x = 20 H, the 
dependence of K and L on these ratios is seen to be much stronger 
(Figure 3.1b). 

Considering first the region directly behind the block where 
the computation is insensitive to the value of C^/Cp and C^/C^Cj^, it 
is reasonable to assume that the ratios C^/Cp and C^/C^Cj^ approach 
unity. Therefore, 
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1 K/K^ 2 1 2 L/Lq 

(b) Profile at 20.0 H behind a rectangular block 

Figure 3.1 Effect of the ratio of constants on the distribution of 
K/Kq and L/Lq behind a rectangular block [11] 




(3.19) 
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(3.20) 


Equations 3.19 and 3.20, of course, are consistent with Equation 3.18. 
Under the above assumptions, the ratio of turbulence production 
rate to turbulence dissipation rate is thus determined solely by the 
local turbulence properties, K and L, and the mean flow deformation; 
see Equations 3.16 and 3.17. 

As the turbulence reestablishes its fully developed state, the 
large eddies will break down through intertial interaction and transfer 
energy to the smaller eddies [28]. As the eddies become smaller, 
energy dissipation due to viscous effects becomes more and more impor- 
tant. Hinze [28] indicates that the turbulence should be characterized 
by the relative rate of change in turbulence energy (DK/Dt)/K. Because 
the rate of change in eddy size DL/Dt is dependent on the local turbu- 
lence properties, this rate of change can be assumed to reflect the 
condition of the local turbulence level, K. A characteristic parameter 
of turbulence can, therefore, be defined as; 


j ^ DK/Dt 
^ ( DL/Dt 


(3.21) 


It is noteworthy that the dimension of appears to be the same as that 
of frequency, (time) . An eddy frequency associated with vorticity 
fluctuation has long been recognized as an important parameter of 
turbulence [3]. 
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Now consider the range where the turbulence is characterized by 


P|^ and The factor can be expressed as 


T 


f 



f C 

V R 



(3.22) 


Based on the assumption for deriving Equations 3.19 and 3.20, the ratio 
of coefficients is again assumed to approach unity, i.e., 

C = 1.0 (3.23) 

y R 

From Equations 3.18, 3.19, 3.20 and 3.23, the relationship between the 
constant coefficients can be written as follows: 



Cp = /r (3.24) 

E y 


It is interesting to note that the value of for smaller eddies where 
and P|^ are dominant becomes 



(3,25) 


which becomes /K/L when Equation 3.24 is introduced. This is consistent 
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with the assumption throughout this section that the relative turbu- 
lence properties; Equations 3.16, 3.17 and 3.22; are insensitive to 
the ratios of constant coefficients in regions of high rate of turbu- 
lence production. Equations 3.24 and 3.25 imply that this, assumption 
is also valid in regions of high rate of turbulence dissipation. 

The value of can be determined from the concept of a constant 
shear layer where the length scale, L, is proportional to the distance 
from the solid wall, i.e., a mixing length hypothesis [11]. 

C = /FM = V*//K (3.26) 

where t is the wall shear stress and is the friction velocity. 

Equation 3.26 indicates that the value of depends on the flow being 
investigated. 

Gosman, et al . [2] propose that the value of can be taken as 
an asymptotic universal constant in high turbulence flows where the 
turbulence Reynolds number = /KL/v is large. The value of 
determined from Equation 3.26 is taken as the asymptotic value when 
applied to the high turbulence level flows considered in this study. 

For atmospheric flows, = 0.416 is proposed [11]. Substituting 
this value into Equation 3.24, the following constants result. 

Cn = C = 0.416 
D y 

= 1.55 

= 0.645 (3.27) 
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These values are comparable to the previous work [11], where 

C = Cn = 0.416, Cn = 1.44 and Cc = 0.60 where used, 
y D K t 

2 

Launder and Spalding [3] indicate that the ratio V* /K lies 

between 0.25 and 0.3 in various flows. A wider range, however, is 

reported by Zeman and Tennekes [35] who give 0.226 to 0.34. Field data 

for atmospheric boundary layers generally show smaller values of 
2 

V* /K as compared with results from laboratory flows. This implies 

that the turbulence generated in the laboratory is smaller than that 

2 

generated in the real world. The ratio V* /K for an atmospheric sur- 
face layer is considered to be between 0.17 [36] and 0.26 [35]. The 

numerical values of the constant coefficients of Equation 3.27 corre- 

2 

spond to the case of V* /K = 0.17. No significant difference of 

2 

numerical solutions has been observed when V* /K = 0.25 is employed. 

Consider the decay processes of turbulence behind grids. 

Launder and Spalding [3] propose that the value of is half of the 
ratio Cq/C^ if the K-L model is employed. With the aid of Equation 
3.19, one can easily show is 0.5, which is comparable to the present 
predictions, Equation 3.27. In order to verify the capability of the 
present model for calculating turbulent flows, atmospheric flow 
over surface objects with different geometries, i.e., forward steps, 
rectangular blocks, fences, and semi elliptical shaped bodies, is 
investigated. Comparison is made with wind tunnel data and full-scale 
field measurements. Good agreement betv/een the analytical solutions 
and the experimental data is obtained and is shown in Section 6.0. 
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4.0 FUNCTIONAL EQUATIONS IN GENERAL ORTHOGONAL 
CURVILINEAR COORDINATES 


The governing equations of motion for the turbulent atmospheric 
boundary layer under neutral stability conditions, Equations 3.1 , 3.2, 
3.9 and 3.14, are formulated in orthogonal curvilinear coordinates in 
this section. Considering a general orthogonal coordinates system; 
u-jsU^jU^j as illustrated in Figure 4.1; a differential line element 
in this system, d^, can be related to the Cartesian coordinates; 
x,y,z; by 

dl = e dx + e dy + e dz 
X y z 

where e^, e^ and e^ are unit vectors in the x, y and z directions, 
respectively. 
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The magnitude of the vector ds can be expressed as 
(ds)^ = (h^du^)^ + (h 2 du 2 )^ + (h^du^)^ 
where the metric coefficients, h.'s, are defined by 


= 
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The gradient is defined by 

„ _ ^1 3 I ^2 8 ,^33 

h-| 3u^ h^ 3U2 h^ 9u^ 

where e-j , 62 and e^ are the unit vectors in the directions 1, 2 and 3, 
respectively. 


4.1 Continuity Equation 


The steady state conservation mass law for incompressible flow 
can be written in vector form as: 

V • ^ = 0 (4.1) 


In the orthogonal curvilinear coordinates system, Equation 4.1 appears 
as 


-> 

V = 




3u^ ^^ 2 ^ 3 '^!^ 9U2 


(h 3 h^V 2 ) 


9u. 




= 0 (4.2) 
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4.2 Conservation of Momentum Equations 


According to the D'Alembert principle, the forces acting on a 
control volume are balanced by an inertia force. The inertia force, 
in turn, represents the inflow of momentum to the control volume 
(Newton Law of Motion). 


Inertia' 

Forces 


+ 


Body ' 
Forces 


■ + 


"Surface" 

Forces 


= 0 


(4.3) 


The terms involving inertia forces and surface forces must be 
examined in the curvilinear coordinate system. In order to simplify 
the analysis, choose the curvilinear coordinates such that the constant 
u^ surfaces are flat planes (Figure 4.2). Therefore, the fluid motion 
in the directions of u-j and makes no contribution to the inertia 
forces in the direction of u^- The center of curvature, 0, for a 
coordinate, described by a position vector r measured from the 
original of the Cartesian coordinates system, lies on a constant u^ 
plane. The equations of flow written in the general orthogonal 
coordinates system are still very complex since the trace of the center, 
0 (Figure 4.2), may be described by a very complex function in the 
x-y-z space. In order to reduce the complexity in the computational 
procedure, the coordinate system used in this study is taken as 
axisymmetric (Figure 4.2). This implies that the u^ coordinate repre- 
sents the angle of revolution, y, about the axis of symmetry from a 
given reference plane, e.g., the x-z plane. Referring to Figure 4.3, 
the inertia component of the rate of momentum change per unit volume, I, 
after neglecting the higher order terms can be expressed as follows: 
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Center of Curvature of ds^ 

Figure 4.3 Velocity vectors of a fluid element shown in a constant 
u^ plane (reference Figure 4.2). 
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where 3 is the angle between the radius of curvature r 2 and the axis of 
symmetry and a = (tt/2) - 3. Substituting 


dt 


and 


dt 


into Equation 4.4, the inertia forces per unit volume, I, become 


~t -> 
I = -pe. 
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Considering the steady state and substituting Equation 4.2 into 
Equation 4.5, the inertia terms in direction 1 are expressed as 


In = 


-P 




9u. 
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pV2V^ pV. 


- + 
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pV. 


sin 3 


(4.6) 


The terms in the first bracket on the right-hand side of Equation 4.6 
can be written as 
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The radii r, and are given by [2] 


1 ^ 1 ^ 


(4.8) 
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(4.9) 




3U. 


Substituting Equation 4.8 into 4.7 and considering the coordinates 
system in which the metric coefficients h-j and h 2 are invariant with 
respect to the direction 3 coordinate. Equation 4.7 becomes 



V- (^h^V^) - 



3h^ 

9U-J 



f- sin e 

3 


(4.10) 


Expressions for the momentum change in directions 2 and 3 are similar 
to that given by Equation 4.10 and are included later in the complete 
equations of motion. 

The surface forces in Equation 4.3 include contributions from 

the fluid pressure and from the shear stresses due both to molecular 

viscosity and to turbulent fluctuations. The contribution of fluid 

pressure to the force per unit volume, ^ , is 

P 

Fp = -VP (4.11) 


The surface stress T can be expressed in general form: 





where 

Tj = ; j = 1,2,3 and i = 1,2,3 
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The component T-. represents the shear stress acting In the ith direc- 
tlon on a surface perpendicular to the jth direction. The components 
making up T.. are expressed as follows [2]: 
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(4.12) 


where p is the effective viscosity defined as 


Pe = V + 


The contributions of the shear stresses to the force on a fluid 

element follow a similar derivation as that for the inertial terms given 

in Equation 4.10. Therefore, the general expression for the momentum 

balance given by Equation 4.3 reduces as follows (see also Reference [2]) 
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Direction 2: 
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Direction 3: 


V . [h3(VV3 - T3)] + X ^ = 0 


(4.15) 


where g.| , 92 and are the components of body force in directions 1 , 
2 and 3, respectively. 

4.3 Transport Equations for Turbulence Kinetic Energy and Length 
Scale 

The governing equations for the kinetic energy, K, and the 
turbulence length scale, L, can be stated as 
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[Convection] = [Diffusion] + [Sources] + [Dissipation] 


In vector notation, these equations are expressed as follows: (see 

also Equations 3.9 and 3.14) 

^ • VK = 1 V • (y^.VK) + P|^ + (4.16) 

V • VL = 1 V • (y^vL) + ''' \ (4.17) 


The expressions and were defined previously in Equations 3.8 
and 3.12, respectively. However, the expressions Pj^ and in the 
curvilinear coordinates system will be discussed more thoroughly below. 
Gosman, et al . [2] derive the form of P|^ and in a curvilinear 
coordinates system from the following consideration: 


Total Shear-Work" 
on Moving Fluid 


Shear-Work Contribution 
to the Kinetic Energy of 
Mean Motion 


Shear-Work Contribution' 
to the Turbulence 
Kinetic Energy 



The total shear-work rate per unit volume can be written as 
”shear = ’ ’ ^ V2 + 


The shear-work rate is made up of W and of W4., which contribute to the 

m t 

kinetic energy of mean motion and turbulence kinetic energy, respec- 
tively. Details of the expressions W^ and W^ are given in Reference [2]. 
Since the present turbulence model is directly concerned with the term 
W^, it is necessary to write out W.^. in the following form [2]: 
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Referring to Equations 3.9 and 3.14, the production rate of the turbu- 
lence kinetic energy per unit mass of the fluid is therefore 


P. = v.W. = C /K L W. 
K t t p t 


and the dissipation rate of L is 




(4,20) 


(4.21) 


Substituting Equations 3.8, 3.12, 4.20 and 4.21 into the corresponding 
terms of Equations 4.16 and 4.17, respectively, the final form of 
Equations 4.16 and 4.17 appears as follows: 


7 • VK = V • (v^VK) + C^/K L W 

V • VL = V • (v.VL) - C Cp — 
t ^ R /K 

where is given in Equation 



+ Cp/K 

4.19. 


(4.22) 

(4.23) 
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4.4 Vorticity and Stream Function Equations 

The procedure for obtaining the velocity components from the 
momentum equations. Equations 4.13 through 4.15, involves the solution 
of the Poisson equation for pressure [37] so that the continuity equa- 
tion, Equation 4.1, is satisfied. For two-dimensional flow situations, 
a stream function, ip, and a vorticity, co, are frequently used to reduce 
the number of variables from three (two velocity components and pres- 
sure) to two (ip and oj) where no prediction of the pressure distribution 
is needed. The vorticity, w, and stream function, ifj, variables are 
used in this study and are formulated as follows. 

Defining the stream function and vorticity with the relationship 


w = J_ -li. 

1 phg BUo 


(4.24) 


^ ^ It 

2 ph-j 3u^ 


(4.25) 


(jO = 


h^h2 


^^ 2^2 ^^ 1^1 


3Ut 3u. 


(4.26) 


the corresponding stream function and vorticity equations are [2]: 


V 4) = -p(jo 

and 


(4.27) 
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For two-dimensional flows, the vorticity vector, w, is always 
perpendicular to the plane of the flow and thus behaves as a scalar. 
Gosman, et al . [2] show that Equation 4.28 can be written as: 


V . (p^co) = V • [V(PgW)] + 


(4.29) 


where S is a source term to be discussed later. Using two vector 

(jO 

operators defined as follows [2]: 


V 


9 . ^ 9 

P — + P 

X 9x z 9z 


V = -e -^ + e ~ 
^ X 9z z 9x 


(4.30) 


where the unconventional symbol V is used in this report for writing 
convenience; and e^ are the unit vectors in the direction of the 
horizontal coordinate, x, and of the vertical coordinate, z, respectively 

The term S is expressed as [2,11]: 

(jO 

= 2[v(6j^ • ^) • ?(e^ • V yg) + V(e^*f) • V(e^ • V y^)] (4.31) 

The expression of Equation of 4.31 in orthogonal curvilinear coor 
di nates is (see Appendix): 
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'ae 

3Ut au« 
L 1 2 

90 
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3^1 ' 
au-j 

* J 


where 


(4.32) 


= x!!!e 
^1 au^ 


^2 




and 0, the angle between e^ and e^ , isdefined in the Appendix, Figure A.l. 

The above equation implies that if the effective viscosity y„ 

s 

is uniformly distributed, i.e., the derivatives of y^ with respect to 
both u-j and U2 are zero or are negligible compared with the other 
terms in the vorticity transport equation. Equation 4.28, the term 
can be wholly omitted from Equation 4.29. 


4.5 Analysis of Flow Parameters 

Since the problem considered here is atmospheric motion over 

bluff bodies, one can choose the characteristic body height, H, and the 

undisturbed upwind velocity at this height, Vn, as the characteristic 

H 

length scale and velocity scale, respecti vely. 

Defining the dimensionless variables 
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H 


V 



" pHV. 


(4.33) 


the resulting dimensionless form of the governing equations (Equations 
4.27, 4.29, 4.22 and 4.23) becomes: 


1 

9 

[^2 


+ 3 



^1^2 

9U.J 

’’l 

9u.| 
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9U2 

h2 2 

V. J 



= -OJ 


(4.34) 
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(4.35) 
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(4.36) 
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h2 2J 


- R.|Fi2 Sj_ = 0 


(4.37) 
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Considering that = (y^ + y)/pHVj^, one can rewrite iJg as 




1 

Re 


(4.38) 


where 


Re 


y 


It is obvious that for turbulent flows where the Re number is very 

~l/2~ 

large, y reduced to C K L and the effect of the molecular viscosity 

c y 

on the flow field can be ignored. The flow is thus insensitive to the 
Reynolds number based on molecular viscosity. 

In a neutrally stable atmospheric boundary layer, the wind 
profile over flat homogeneous terrain may be described by a logarithm 
law [4]: 


V = 


V* 

— In 

K 



where 




(4.39) 


The parameters V* and Zq are the friction velocity and surface roughness 
length scale, respectively, and von Karman's constant k is taken as 0.4. 
The friction velocity V* appearing in Equation 4.39 can be thought of 
as a measure of the turbulence level in the approaching flow (see 
Equation 3.26). 
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Examining Equations 4.35 through 4.39 and being cognizant of the 

above argument, the basic parameters involved in low-level atmospheric 

motion, where Coriolis effects are negligible, are found to be p^, V*, 

~l/2~ 

and Zq, where = C^K ' L. The first is a measure of turbulence trans- 
portation due to turbulence fluctuations, whereas the second and third 
characterize the turbulence level in the undisturbed upstream flow. 
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5.0 NUMERICAL ALGORITHM 


The governing equations. Equations 4.34 through 4.37, given 
in Section 4.0, are solved by utilizing the numerical procedure 
developed by Gosman, et al . [2]. An outline of the procedure is 
given in Section 5.1. Theoretically, the procedure can be applied 
to any curvilinear coordinate system, however, a semi elliptical 
cylinder on a plane surface, discussed in Section 5.2, is chosen in 
this study. This geometry simulates a hill in natural terrain for 
which parametric variations of the aspect ratio can be carried out. 

The boundary conditions for the dependent variables w, $, K and L 
are given in Section 5.3. In Section 5.4 the accuracy of the numerical 
results IS discussed. 

5.1 Derivation of Finite Difference Equations 


Following the method of Gosman, et al . [2], Equations 4.34 
through 4.37 can be written in the general form: 


- 


■ 


' ho 


f h 1 


3uJ 

■ 3uJ 

- 

9 

9u-| 

' h 

1 "■] ' J 

3U2 

h2 ^“2 J 



advection ^ ^ diffusion 


+ h-j ii2 d = 0 (5.1) 

■-source -■ 

where 4 > stands for w, K or L. ' The corresponding functional coeffi- 
cients a, b, c and d are listed in Table 5.1. 
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Table 5.1 Coefficient Functions of Equation 5.1 


<p 

a 

b 

c 

d 

$ 

0 

1 

1 

-ai 

to 

1 

1 


-S 




^e 

0) 

K 

1 


1 

“^K 

L 

1 


1 



The general differential equation. Equation 5.1, is replaced by a finite- 
difference equation in the numerical computation procedure. The finite- 
difference equation is obtained by integrating Equation 5.1 over a 
typical cell appropriate to a central node p surrounded by nodes E, W, N 
and S, Figure 5.1. The resulting finite-difference equation is [2,38]: 
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L — source > 


where 
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S = (“l)e ■ (“l)w = ^“2 = - (“ 2)5 

= (®l>e - •> (^®2^ns = (®2^n ‘ 


s^ and $2 are new independent variables which are connected to the 
u^-U 2 coordinates by the relations: 

ds^ = h^du^ ; dS2 = h2dU2 


The approximation for the advection terms in Equation 5.2 employs an 
"upwind differencing" technique. This is a one-sided rather than a 
centered-space differencing, where the scheme is backward differencing 
when the velocity is positive and forward differencing when it is 
negative. This formulation of the first-order terms gives greater 
numerical stability than can be obtained with central differences [32]. 
Using the upwind differencing method, the first advection terms is 
approximated by [2]: 



where aJ = 

^ne ^ l^^NE 

Thus the value of (j^g in Equation 5.2 appears as follows: 
4>g = <\>p if A$ > 0 
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■'■f' < 0 

The term (9i|;/9u2)g has been approximated by 


' dip ' 
9Uo 

I 




The remaining advection terms are obtained similarly and the resulting 
finite form of the advection terms is: 


[Advection] ~ 


P^P 


I A-'f-i 

i=NSEW ’ ^ 


(5.4) 


where 


i=NSEW ^ ^ ^ ^ ^ ^ 


S^^SEN I^SEnI^ 

\ I^NWsI^ 

" S^^ENW I^ENwl^ 

" 8*^^WSE I^WSeI^ 

Ap = Ae + A^ + An + A3 


+ A^cD^ + A^4>^ 

’^SEN " %E " ^NE 

^NWS " ^NW ' ^SW 

^ENW " ^NE " 

^WSE " %W ■ ^SE 




The above definition of the symbol I’s used throughout this study. 

Roache [37] indicates that the application of "upwind differencing" 
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technique may introduce an artificial viscosity. This pseudo viscosity, 
however, can be reduced by decreasing the grid size used in the numerical 
programming. The influence of the pseudo viscosity on the flow field is 
discussed in Section 5.3. 

A central -difference method is used to approximate the diffusion 
terms in Equation 5,2. It is assumed that the b's and c's can be 
approximated with a linear variation, i.e.. 


3 , , , 

~ 7 ^ 7 

^^1 ®1.E ■ ®1,P 


(5.6) 


The (Asi)g^ and (AS 2 )p 5 terms appearing in Equation 5.2 are approximated 
as 


(Asi ) 


ew 




(5.7) 


and 


(As 2 ) 


ns 


2^^2,N ■ ^2,S^ 


(5.8) 


respectively. Writing the remaining b's and c's in the forms of Equations 
5.5 through 5.8, the resulting finite-difference form becomes: 


[Diffusion] = T 

i=wkw 


B^(c()>)^ 


-Bp(c*)p 


(5.9) 


where 
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_ “’e '’p ®2;,N ■ ^2,S 

E - 4 ~ 

®1(E " 5l.P 


^ ^ ‘’P ^2,N ~ ^2,S 

®1.P ■ ®1.W 

_ + bp S-,^£ - 

’ s - ; 

^2,N ®2,P 

^ ‘’S '*' *^P ^1,E ' ^1,W 
^2,P ■ ^2,S 

= Be + By + Bn + Bs 


Introducing Equations 5.7 and 5.8 into the source terms, the following 
finite-difference equation is obtained. 


[Source] ^ ‘^p^^l.E ' ^1,W^^®2,N ‘ ®2,S^ 


(5.10) 


Summarizing, the finite-difference form of Equations 5.4, 5.9 and 5.10 
for integration of Equation 5.1 is 


“ I 

^ i=NSEW ^ ^ 


+ D 


(5.11) 


where 


C. = (A. + B.c^.)/™ 
D = 'dpVp/zAB 
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’ ^2,S^ 


IAB = I A| + r B. 
i=NSEW ^ P ^ 


Because of the nonlinear characteristics of the resulting finite- 
difference equation, Equation 5.2 is solved by an iterative 
successive substitution technique. 

5.2 Elliptical Cylinder Coordinate System 


The orthogonal- curvil inear coordinates system used in the study 
is an elliptical cylinder system. This choice allows one of the 
coordinate lines to represent the solid boundary of a semiel 1 iptical 
body which in turn is assumed to simulate a hill. 

Mathematically, a two-dimensional semielliptical surface, A-B 
as shown in Figure 5.2(a), can be defined by the following equation: 

2 2 

% + % = 1.0 (5.12) 


where D is one-half the longitudinal length of the semiel 1 iptical body 
and H is the body height. In two-dimensional elliptical coordinates 
(u-|, U 2 ), Figure 5.2(b), the coordinates x and z can be expressed as 

X = -a cosh U 2 cos u-j 

z = a sinh U 2 - sin u-j (5.13) 

where a is one-half the length between the two focus points of the body; 
U 2 ^ 0 and 0 < U-J < 2 tt. 
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(a) Schematic diagram of the numerical coordinates used 



(b) Two-dimensional semielliptic coordinate system 

Figure 5.2 The numerical and physical coordinate system for calculating 
the problem of flow over semiellipses 
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From Equations 5.12 and 5.13, it is evident that 


H = a sinh u. 


D = a cosh u. 


(5.14) 


and therefore, a value [u 23 |^ representing the solid ellipse boundary is 
determined by the following equation: 


" coth’^ ^ 


(5.15) 


From Equations 5.14 and 5.15, the dimensionless form of Equation 5.13 is 


X = -cosh cos u-|/sinh (u2)j^ 
z = sinh u^ sin u-j/sinh ( 02 )]^ 


(5.16) 


The trace of the curvilinear coordinate surfaces of constant u-j and U 2 » 
respectively, on the x-z plane shown in Figure 5.2(b) are obtained from 
Equation 5.13. They are confocal ellipses and hyperbolas, respectively. 
The dimensionless metric coefficients h-j and h^ are as follows: 


h^ = a/sinh2 U 2 + sin2 u^ 


h2 = 


where (see Equation 5.14): 


a = 1/sinh [u2]|^ 
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5.3 Numerical Coordinate System 

The origin of the numerical coordinates is situated at the up- 
stream corner of the semi elliptical body, node A in Figure 5.2(a) with I 
indexing in the positive u-j direction and J indexing in the positive 
direction. In the iteration process, the field was swept from left to 
right beginning at the outer boundary and proceeding in the decreasing 
J-direction. A total of 2,501 grid points (I x J = 61 x 41) are con- 
structed in the flow field. The variable mesh used is determined by 
the following equation: 

[u^lj = (I - 1 )tt/60 ; I = 1, 2, 3, ••• 61 

[U2]j = [U2]j=i + I ; J = 2, 3, 4, ••• 41 (5.17) 

n =2 

where 

[U 2 ]j=l = coth'’ p 

and C-j and are constants, A value assigned for C-j , then, represents 
the distance of the first grid point to the obstacle through Equation 
5.16. A value greater than unity assigned for 62 * however, is used to 
stretch the grid size and therefore changes the distance of the outer 
boundary from the obstacle. For D/H = 2.0, Equations 5.16 and 5.17 show 
that the outer boundary is a distance of approximately 30 H from the 
obstacle if = 0.025 and C 2 = 1.05. The distance becomes 45 H if C 2 

is changed from 1.05 to 1.055. The distance from the first grid point 
to the obstacle is approximately 0.05 H in these cases. A distribution 
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of grid points in the x-z plane is shown in Figure 5.3 for the case 
D/H = 2.0, = 0.025 and = 1.05. 

5.4 Boundary Conditions 

The dimensionless form of the boundary conditions is discussed 
in this section. The outer boundary conditions are determined from the 
undisturbed flow conditions, whereas the obstacle and the wall boundaries 
are determined from the local equilibrium conditions. 

5.4.1 Outer boundary . It is assumed in the present study 
that the outer boundary. Figure 5.3, is located far away from the semi- 
elliptical body on the half plane so that the flow conditions at the 
outer boundary are not disturbed by the presence of the obstacle in 
the flow field. For atmospheric flow as considered in the study, the 
undisturbed wind velocity profile is the logarithmic velocity profile 
given by Equation 4.39. It is (Equation 4.39): 



K 



where von Karman's constant k is taken as 0.4. In the Cartesian co- 
ordinates system, the undisturbed stream function $(z) is obtained by 
integrating this velocity profile over the height of the solution region. 


^{z) = 


V^z. 


V dz = 


(1 + z^) £nO + z^) - z^ 


(5.18) 


59 



Figure 5.3 Distribution of grid sizes for D/H = 2.0,. C. = 0.025, 
=1.05 ' 


where 1s the ratio z/Zq* For a grid node, Q, on the outer boundary. 
Figure 5.3, the height, Zq, can be determined from its coordinates 
[(u-j)q, (o 2 )q] by the following equations (see Equations 5.13 through 
5.15): 

_ s1nh(u2)Q sin(ui)Q 
^ sinh[coth"^ (D/H)] 


where (u])g and (u 2 )q are the values of u-j and U 2 at the position of 
the node Q, respectively. At the node Q the boundary condition of the 
stream function T/jg is therefore (see Equations 5.18 and 5.19): 


" ~1T 


1 


In 


1 + ^ 


!a 



(5.20) 


Similarly, all other boundary conditions of 41 for the nodes on 
the outer boundary are obtained. This procedure is also used to obtain 
the outer boundary- conditions for w, k and L, respectively. Thus, only 
the undisturbed conditions for w, ic, and L, respectively, of the 
atmospheric boundary layer written in terms of the z are presented in 
the following. Equation 5.19 is used to convert (u-|,U 2 ) to the (x,z) 
coordinate system. 

In the undisturbed atmospheric boundary layer, the dimensionless 
vorticity ai is written as follows: 


w(z) = - 


3V 

8z 


V 




K Z 




(5.21) 
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The boundary condition for the turbulent length scale is based; , 
on a mixing-length hypothesis. In the neighborhood of the wall, the 
mixing length is generally assumed to be linearly related to the normal 
distance from the surface. For a very rough surface, the mixing length 
hypothesis does not go to zero but approaches a size on the order of the 
roughness length, z^. The following relationship is assumed; 

L = k(z + z^) (5.22) 

The boundary condition for the turbulence kinetic energy, K, is derived 
from the assumption of constant shear layer. In terms of the Prandtl- 
Kolmogorov formula. Equation 3.5, the shear stress, t, of the undisturbed 
atmosphere is written as follows: 

T = pvj = c^p/i< l[|| 

Substituting Equations 4.39 and 5.22 into the above equation, the 
following equation is obtained: 

'V J 2 

K = -^ (5.23) 

. 

5.4.2 Wall boundary . The wall boundary is assumed to be a 
solid wall; therefore, no flow passes through it. Thus, the wall boun- 
dary forms a streamline which is assigned the value zero. 

$ = 0 (5.24) 
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Applying the nonslip conditions at the wall; i.e., V-j = 0; V2 = 0; the 
w at the wall becomes (Equations 4.26 and 4.33): 


(jO 


wall 


3V, 

i 

3u, 


u 


'wall 


In the present study, the wall boundary is constructed by u^ =0, 
u^ = 7T and u^ = (^2)15 (see Figure 5.2 and Equation 5.15). Considering 
that 9V^/3U2 = 0 along u^ = 0 and u^ - ir, and that 3V2/3u-j = 0 along 
^2 ~ ^^2^b* dimensionless vorticity along the wall can be written as 
follows: 


03 


wall 


3V. 


3u- 


wal 1 


if u^ = 0 or = TT 


and 


3V. 


0 ) 


wal 1 


3u. 


if u« = (u^) 


2'b 


wal 1 


Since the above equations can be written in a general form when the 
velocities V-j and V2 are replaced by Equations 4.24 and 4.25, respecti vely , 
only the expression of wall vorticity along U2 = (u2)|j» terms of 

stream function (see Equation 4.24) is discussed in the following: 


“w 


h2 




(5.25) 
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To express in terms of known variables, the stream function around 
a nearby wall point, point NW in Figure 5.2a, is expanded in a Taylor 
series : 


^NW 3u2 


AU„ +1^ 


7 7 

^ 802^^ 


+ H.O.T. 


w 


(5.26) 


The second term on the right-hand side of Equation 5. '26 is zero due to 
the nonslip condition. From Equations 5.25 and 5.26, the vorticity at 
node W on the wall boundary; can be approximated by the following 
equation if the higher order terms of Equation 5.26 are neglected. 


h2(AU2)^ 


(5.27) 


Assuming that the concept of a constant shear layer is applicable in the 
region very close to the wall, the kinetic energy of turbulence, K, and its 
length scale, L, take the following form, respectively [11]. 



I Vi j 

L = K (5.29) 

5.5 Accuracy of Numerical Solutions 

The accuracy of numerical solutions is generally affected by the 
imposed outer boundary condition^and the grid sizes used in the solution. 
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Theoretically, the limit of the outer boundary is at infinity. In a 

practical analysis, however, this cannot be imposed; therefore, in this 

study the conditions at infinity are applied at a large distance, 

from the elliptical body. The value of z is determined in such a way 

rriaX 

that further increase in its value does not significantly effect the 

numerical solutions. The value of is determined by numerical experi- 

max 

ment. In the case D/H = 2.0, solutions obtained from z , = 30 H and 

max 

45 H, respectively, are compared in Figures 5.4 through 5.7. This com- 
parison indicates that the difference between these two solutions is 
small (about 5% of the local value). The solutions for the case D/H = 2.0 
shown in Section 6.0 are therefore computed with z^^^ = 45 H to assure 
good compliance with the boundary conditions applied at the outer boun- 
dary, Equations 5.18 and 5.21 through 5.23. However, a value of 

= 30 H would provide the necessary accuracy if desired. In the 

niaX 

present study, the numerical computation for D/H = 1.02 is carried out 

for z = 30 H. 
max 

In addition to the outer boundary condition, the grid size used 
in a numerical procedure can significantly affect the accuracy of the 
solutions. Two cases are used to investigate the effect of the distri- 
bution of grid size on the numerical solutions. The first case is 
C-j = 0.05 and C 2 = 1.02 in Equation 5.17, and the second case is 
C-j =0.05 and C 2 = 1.05. As mentioned in Section 5.3, the number 
of grid points constructed for numerical computation is 2,501 points 
(I X j = 61 X 41), which yields z^,^ = 30 for both cases mentioned 
above. The number of grid points in the wall region of the second case 
is approximately twice the number of the first case. In the region 
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= 45 
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max 


z 

= 30 
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max 



Figure 5.5 Distribution of turbulence kinetic energy, K, for D/H = 2.0 









near the outer boundary, more grid points are constructed for the first 
case than are constructed for the second case. Again, the effect of 
the grid size on the numerical solutions is small for the cases compared. 
Figure 5.8. Based on the above discussion, it was assumed that the grid 
distribution given by Equation 5.17 with C-j = 0.025 and C 2 = 1.055 
would provide meaningful results. Therefore, this distribution is used 
throughout the study. The application of C-j = 0.025 and C 2 = 1.055 
reduces the grid size to 0.05 H at the wall region and approximately 
7 H at the region near the outer boundary if D/H = 2.0. For D/H = 1.02, 
the grid size at the wall region and at the outer boundary region is 
0.025 H and 5 H, respectively. The distance from the outer boundary to 
the obstacle is approximately 45 H if D/H = 2.0 and is approximately 
30 H if D/H = 1.02. 

As mentioned in Section 5.1, the "upwind differencing" technique 
used in this study causes a pseudo viscosity, the value of which depends 
on the grid size [37]. Since the grid size near the obstacle is very 
small (not greater than 0.05 H), the effect of the pseudo viscosity is 
considered to be small in that region. An analytical method of 
determining the error caused by the pseudo viscosity effect is not 
available for nonlinear equations; the reader must therefore make his 
own judgement of this effect from the grid distribution shown in 
Figure 5.3. From previous experience, the author believes that this 
effect will not cause significant inaccuracies in the results. 

The convergence criteria for the numerical computation is 




n-1 

P 


< 10 


-4 


(5.30) 


max 
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Figure 5.8 Velocity profiles on top of a semiellipse for 
D/H = 2.0 and = 0.005 
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where {j> is the dependent variable at grid point P, and the superscript n 
denotes the nth iteration. Equation 5.30 insures good convergence in 
areas of small <})p values. The present study was run on an IBM 370/3031 
computer. Execution of the present computer code requires approximately 
256 K bytes of main storage. Each iteration of the code requires approxi- 
mately 9.5 sec on the central processing unit of the IBM 370/3031. The 
computational time required for convergence of the cases run in this 
study is approximately 180 minutes for each case if the undisturbed flow 
conditions, i.e., outer boundary conditions described in Section 5.4 
are the initial conditions. The computation time is reduced to 120 
minutes if the converged solutions of a case are used as the initial 
conditions of other cases. 


72 



6.0 RESULTS AND DISCUSSION 


The numerical results obtained in this study represent solutions 
of atmospheric flow over a two-dimensional semiel 1 iptical body. The 
effects of aspect ratio, of distributions of surface roughness, and of 
molecular Reynolds number are discussed in this section. The numerical 
results are compared with the available experimental data. Since most 
currently available experimental data are for flow over rectangular- 
shaped bodies, i.e., fences, blocks and forward-facing steps, numerical 
results computed with the present model in a Cartesian coordinates system 
over these types of geometries are first tested. The agreement between 
experiment and theory is excellent. Since the computation algorithm used 
to compute flow over the semiellipse is the same as that for Cartesian 
coordinates, the present computed results for flow over semiellipses are 
believed to be good even though there is insufficient data available 
for a valid comparison. The flow separation created by a rectangular 
geometry, however, is quite different from separation caused by a 
semi ell ipse. The differences in the mechanism of flow separation are 
discussed in Section 6.5. 

6.1 Characteristics of Flow About a Two-Dimensional Semielliptical 
Obstacle 

The results discussed in this section are obtained for a semi- 
elliptical body with aspect ratio D/H = 1.02. The outer boundary is 
situated approximately 30 H away from the body. The dimensionless sur- 
face roughness scale, Zq = 0.005, is uniformly distributed along the 
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wall and over the semi ellipse. Substituting Zq = 0.005 into Equation 
4.39, the dimensionless friction velocity is obtained as V* = 0.075. 
From Equations 3.27 and 5.23, the undisturbed upstream turbulence 
energy is therefore K = 0.0325. The molecular Reynolds number. Re, is 
assumed to be large enough so that its effect on the solutions is 
negligible, as discussed in Section 4.5. This assumption is reason- 

5 

able when Re is greater than 10 , as found from the present numerical 
study. The effect of Re is discussed later. 

The streamline pattern for Zq = 0.005 and D/H = 1.02 is shown 
in Figure 6.1. The flow is seen to converge as it passes over the 
obstacles. The convergence of the flow causes the wind to gain speed 
in the region above the obstacle. Figure 6.2. The constant velocity 
contours in Figure 6.2 show that the local maximum in velocity occurs 
at 0.1 H above the crest of the body and the local minimum occurs at 
1.0 H. Similar phenomena have been observed in full-scale field 
measurements [9]. More discussion of this phenomenon is given when 
the numerical results are compared with the experimental data in 
Section 6.4. 

Figure 6.2 shows very strong gradients occurring in the area 
close to the surface of the obstacle. This implies that a high produc- 
tion rate of turbulence kinetic energy, K, occurs there, as discussed 
in Section 3.0. The production rate of K in that region is stronger 
than its dissipation rate and, therefore, relatively high turbulence 
kinetic energy occurs in this region, as shown in Figure 6.3. Down- 
stream of the obstacle, the region of high turbulence expands grad- 
ually, Figure 6.3, due to the effects of convection and diffusion. A 
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relatively low turbulence area is observed in the region approximately 
1 H to 2 H above the crest, as is expected due to the nearly uniform 
velocity profile. Figure 6.2. In the vicinity of both the upstream 
and the downstream corners of the obstacle, the turbulence kinetic 
energy is very small due to the motion of the fluid being restricted 
by the wall . 

The distribution of the turbulence length scale is shown in 
Figure 6.4. This figure indicates that the length scale L decreases as 
the flow approaches the obstacle. Then L increases rapidly as the flow 
passes over the obstacle. Downstream of the obstacle, the L decreases 
again with the exception that L increases in a region immediately 
behind the obstacle. The rapid increase of Z on top of the obstacle 
implies that the length scale cannot be prescribed in that region by 
a mixing length hypothesis. At 0.5 H above the crest, the length scale 
L is equal to 0.4, which is the same value as L upwind at 1 H above the 
plane surface. The present results indicate that L would be under- 
predicted by a mixing length hypothesis in the region above the 
obstacle. 

Constant vorticity contours are shown in Figure 6.5. Very 
strong vorticity is generated on the surface of the obstacle due to the 
nonslip condition at the wall. The distribution of vorticity in the 
flow field is determined by the velocity gradients. Convection and 
diffusion of vorticity are clearly shown in Figure 6.5. 
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Figure 6.4 Distribution of turbulence length scale, L, over a semiellipse 

for D/H = 1.02 and z = 0.005 
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6.2 Effect of Aspect Ratio 


The effect of variation in aspect ratio on the flow field for 
a nondimensional surface roughness Zq = 0.005 is shown in Figure 6.6. 

The V/V^ isotachs for the aspect ratio D/H =1.02 and 2.0 are presented. 
One observes that the velocity above the crest of the obstacle increases 
with a decrease in aspect ratio. This effect is very strong in the 
lowest 1 H layer above the crest and is felt even up to 10 H, Figure 
6.7. Reference [19] explains that a stronger favorable pressure gradi- 
ent occurs near the leading edge of a smaller aspect ratio semi- 
elliptical body, thus causing higher wind speeds on top of the body. 

The constant contour lines of the dependent variables ifi, w, 

K and L for D/H = 2.0 were shown in Figures 5.4 through 5.7. Comparing 
these figures with Figures 6.1 and 6.3 through 6.5, one observes that 
the vorticity, w, produced near the surface is significantly affected 
by the aspect ratio. However, the effect on the turbulence kinetic 
energy is small. The reason for this is that by definition the vorticity 
is a direct function of the velocity gradient; however, the local tur- 
bulence kinetic energy, K, is a balance between the diffusion, the pro- 
duction, the dissipation and the convection of K. Therefore, the local 
turbulence kinetic energy is dependent not only on the velocity gradient 
but also on the local turbulence level and the surface roughness. The 
effect of aspect ratio on K is small for the cases compared in the present 
study because the surface roughness has been held constant (2 q = 0.005) 
along the wall. Moreover, the undisturbed turbulence kinetic energy also 
remains the same, i.e., K = 0.0325, for the same surface roughness value. 
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6.3 Effect of Surface Roughness 


The effect of surface roughness on the flow field is described 
in two parts: (1) the effect of the local value of assigned to the 
obstacle and (2) the effect of the value of Zq assigned to the upstream 
plane. 

The effect of the local Zq is studied for the aspect ratio 
D/H = 2.0. The value of Zq used to determine the undisturbed friction 
velocity, V*, Equation 5.14, is 0.005 in this case. The value of V* is 

0. 075, and therefore the undisturbed upstream K is 0.0325, Equation 5.23. 
Numerical results are obtained for z^^ = 0.01 and for z^ = 0.005, where 

is the surface roughness scale of the obstacle. For increasing Zj^, 
the velocity decreases in the lower layer directly on top of the 
obstacle. Figure 6.8. The thickness of this layer is approximately 
1 H. The increase of surface roughness causes a larger drag on the 
flow field. This friction drag reduces the wind velocity in the lower 
layer. Figure 6.9, and thus displaces the streamline to a greater 
height. Figure 6.10. Since the kinetic energy lost in the mean wind 
field due to the effect of friction drag is converted into the turbulent 
fluctuations, the turbulence kinetic energy is increased with increasing 
z^, Fi gure 6.11. 

The flow field about an obstacle is influenced not only by the 
roughness scale on the obstacle but also by its distribution along the 
wall. If the roughness scale is increased from Zq = 0.005 to z^ = 0.01, 

1. e., V* increases from 0.075 to 0.087, the turbulence kinetic energy, 

K, in the upstream undisturbed flow increases from 0.035 to 0.044. 
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Figure 6.8 Velocity profiles on top of a semiellipse with different 
surface roughnesses for D/H = 2.0 








The K in the flow field increases correspondingly, as one observes 
in Figures 6.11 and 6.12. From Equation 5.18, the maximum value of 
ijj assigned at the outer boundary of the computed flow region, i.e., 
the value of at z = 45, is 68.4 if Zq = 0.005 and V* = 0.075, and 
is 76.2 if Zq = 0.01 and V* = 0.087. Therefore, more mass flow is 
brought into the flow region computed for the case of larger Zq and 
thus larger velocity in the flow field is obtained. Figure 6.13. 

The velocity profiles on top of the obstacles, shown in 
Figure 6.8, indicate a slower wind speed in the layer 0.2 H above the 
crest for = 0.01 than for Zq = 0.005. This is obviously a result 
of the larger drag force created by the rougher surface. In the layer 
between 0.2 H and 1.0 H, V increases for larger Zq due to t'le enhanced 
turbulence mixing process which carries the upper layer fluid elements 
with high momentum into the lower layer. This interaction augments 
the momentum in the lower layer. 

Figure 6.14 illustrates that the velocity speed-up ratio, 
defined as the local speed at a height z above the crest of the obstacle 
divided by the undisturbed speed at the same height above the upstream 
surface, i.e., S = V(z)/Vq(z), increases with increasing Zq. This 
agrees with the results of previous studies [32]. As explained pre- 
viously in Section 6.2, the aspect ratio has a significant influence 
on the speed-up ratio. From the limited cases discussed, the effect 
of the aspect ratio on the speed-up ratio can be seen to exceed the 
effect caused by the change of surface roughness. Figure 6.14. 


89 







x/H 

: - 0.005 and z. = 0.01 

0 b 

: z^ = 0.01 and z^^ = 0.01 


Figure 6.13 V/V^ isotachs over a semiellipse for D/H = 2.0 






6.4 Verification of Results 


Comparisons of the numerical results with experimental field 
data and with atmospheric boundary layer wind tunnel data are given 
in this section. Since the flow field is significantly affected by 
the geometry of the obstacle, the surface roughness scale and the 
upstream turbulence level (discussed in Sections 6.1 through 6.3), 
only experimental data obtained under well-defined conditions are 
used for comparison. 

Most data available for comparison are for flow over rectangular 
bluff bodies, i.e., steps, blocks and fences. Limited data are avail- 
able for flow over other geometries. In this section, therefore, 
comparison is first made for the case of flow over rectangular bluff 
bodies. For flow over semiell iptical bodies, the numerical solutions 
are compared with the analytical results of Frost, et al . [19]. The 
present solutions are also compared with the results of flow over 
obstacles with geometries similar to semielliptical bodies. These 
cases are the analytical results of flow over a two-dimensional Gaus- 
sian hill [39], the wind tunnel data of flow over two-dimensional 
half- and full-sinusoidal bodies [7], the flow over a three-dimensional 
Gaussian hill [7], and the data of Bradley [9] obtained from field 
studies of flow over a full-scale hill. 

Numerical results for flow over a 90° escarpment computed with 
the current computational model in Cartesian coordinates [32] and the 
experimental wind tunnel data [6] are shown in Figure 6.15. The 
simulated undisturbed upstream flow in the wind tunnel corresponds to 





atmospheric boundary flow over ground with 2 q = 0.1 m if a model scale 
1 :‘300 is assumed [6]. Since the height of the model in the wind 
tunnel study is 50 mm, at 1 : 300 this corresponds to a full-size scale 
body of 15 m. Accordingly, the dimensionless surface roughness 2 q is 
0.007. The numerical results shown in Figure 6.15, therefore, are 
computed with Zg = 0.007 uniformly distributed along the wall. 

The numerical solutions obtained by using the coefficients C^, Cg, 
and given in Reference [32] are also shown in Figure 6.15. The 
results are relatively insensitive to the coefficients used. 

Some differences between the numerical results and the 
experimental data [6] are shown in Figure 6.15. On top of the front 
upper corner of the step, the velocity speed-up ratio, V(z)/Vg(z), is 
slightly over-predicted. In other regions of the flow field computed, 
the numerical results are under-predicted in the region near the 
ground. These observations suggest that the surface roughness scale, 
Zg, used in current computation is higher than the experimental value, 
as stated in Section 6.3. 

Figure 6.16 shows a comparison between the computed value of 
/K/V^ and the experimental data for [6], where is the upwind 

velocity at a height 10 H above the ground and is the r.m.s. value 
of the turbulence fluctuation of the horizontal velocity components. 
Over level homogeneous terrains, many micrometeorological experiments 
suggest that 

ctv 2 “ 0.5 Oy^ 

Ow = 0.8 Ow (6.1 ) 

3 ''l 
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where are the r.m.s. values of the turbulence fluctuations 

of the vertical and the lateral components, respectively [9]. By the 
definition of K, 

^ ~ 2 *^ v -| ^^2 ^ (^* 2 ) 

and the relationships of Equation 6.1, the following approximate result 
is obtained: 

Ov^//K=1.03 (6.3) 

for atmospheric flow over uniform terrain at a neutral stability con- 
dition. However, in actual fact the ratio a^^/ZK is affected by the 
terrain over which the flow passes [40]. As flow passes over obstacles, 
the ratio increases slightly, while the ratio ^v^/'^v-j i^^niains 

approximately equal to upstream value [7,9]. Limited data [9] shows 
that on top of a hill, the ratio increases from its upstream 

value 0.5 to a local value of 0.8. This reduces the value of the ratio 
from 1.03, Equation 6.3, to 0.93. Since the value of the ratio 
Ovi/ZK depends on the terrain over which the flow passes, the comparison 
shown in Figure 6.16 should be made qualitatively rather than quanti- 
tatively. Figure 6.16 shows that the present results of »^/V^ agree 
well qual i tati vely with the experimental data of [6]. The 

analytical results of Zi</V^ obtained by using the coefficients given 
in Reference [32] result in an over-prediction of ZK/V^ of approximately 
five times the present results, in the region next to the wall. It 
is concluded that the coefficients developed in the present study result 
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in an equivalent prediction of the mean wind field and a better pre- 
diction of the kinetic energy field, as compared with those developed 
in Reference [32]. 

Results from the present numerical model applied to a block 
geometry, when compared with experimental wind tunnel data and full- 
scale field data, show excellent agreement in both the mean velocity 
field and the turbulence kinetic energy field [11,26]. It is noteworthy 
that the upstream turbulence level causes significant effect on the 
flow field computed. An increase in upwind turbulence reduces the dis- 
tance, between the block and the flow reattachment point. Figure 
6.17. Reference [41] reports that reattachment of flow over a solid 
fence occurs at a distance 10 H to 15 H measured downstream from the 
fence. This is accurately predicted by the present numerical model in 
the range of V*/V|^ = 0.053 to 0.015, Figure 6.17. 

Figure 6.18 shows typical V/Vq isotachs computed and measured 
for a solid fence. Excellent agreement between the computational 
results and the experimental data [41] is observed. In the recirculat- 
ing flow region, the present results show a region of negative velocity, 
as one would expect. This region of negative velocity, however, is not 
shown in the experimental data. The measurement of velocity using hot 
wire [41] in the region close to the fence, where the velocity is low 
speed and highly fluctuating, can cause error in the measurements. 

Bradley [9] measured wind speed profiles on the top of a hill. 
The terrain surrounding the hill is roughly a uniform surface with 
Zq = 0.005. The hill, however, is covered by a forest. The shape of 
Bradley's hill, along with the shapes of other surface geometries for 
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Figure 6.17 Effect of the turbulence of the approaching wind on the size 
of the wake zone behind rectangular blocks 
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Figure 6.18 V/Vq isotachs downstream of a solid fence, 






which either empirical or analytical data are available for comparison, 
is shown in Figure 6.19. Measured and computed vertical profiles of 
horizontal velocity on the top of the respective hill geometries shown 
in Figure 6.19 are presented in Figure 6.20. The velocity scale Vp, 
used in Figure 6.20 is the upwind velocity at = 0.2 H, where is 
a reference height. 

Great differences between the various analytical models is 
clearly seen from Figure 6.20. These differences are due to the 
different assumptions used in the analyses. The length scale used 
in References [19] and [39] has been prescribed by the mixing-length 
hypothesis. This approach can cause significant differences in the 
analytical results (see Section 6.1). The boundary layer approxima- 
tion of Frost, et al . [19] and the different closure form used by 
Taylor [39] for the shear stress terms of the governing momentum 
equations can also create disagreement between the analytical models. 

The data of Bradley [9] show very low values of V/Vj^ in the 
region z < Zp. This is believed to be caused by the existence of a 
forest on the hill. The existence of the forest increases the effec- 
tive surface roughness and therefore decreases the velocity in the 
region next to the ground (see Section 6.3). An effort was carried 
out to simulate the effect of a forest on the flow field by increas- 
ing the surface roughness scale, Zj^, on the semiellipse. The value 
0.1 was assigned to z^, while the ground surface roughness, Zq, 
remained at 0.005. No converged numerical solutions are obtained 
due to the large magnitude of the abrupt change in surface roughness. 
The numerical solutions for Zj^ = 0.01 and Zq = 0.005 show only slight 
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Figure 6.20 Vertical velocity profiles on the top of the respective 
obstacles shown in Figure 6.19 
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reduction of the ratio V/Vp in the region z < z^. An accurate method 
of simulation of the effect of a forest has not yet been developed. 

Figure 6.20 shows that the analytical results of Frost, et al . 

[19] agree well with the data of Bradley [9]. It is noteworthy, how- 
ever, that the analytical results are obtained by a two-dimensional 
analysis with uniform Zq [19], while the experimental data are obtained 
for flow over a hill with a forest on its top [9]. The wind tunnel data 
for a three-dimensional Gaussian hill [7] are also consistent with the 
computed results of Frost, et al . [19] in the region z > Zj^, Figure 
6.20. As shown in Figure 6.20, the present computed results agree 
excellently with the wind tunnel data [7] of flow over two-dimensional 
obstacles with geometries similar to a semiellipse. In the region 
z < Zp, Taylor's predictions [39] agree with the present computed 
results and then shift to agree with the computed results of Frost, 
et al . [19] if z > z^. Figure 6.20. However, the value of Zq used by 
Taylor [39] is 0.0005 rather than 0.005 used by Frost, et al . [19] and 
by the present computation. 

Although the limited data compared in Figure 6.20 are obtained 
from flow over obstacles with similar geometries, one cannot form any 
conclusions from the agreement (or disagreement) between the analytical 
results and the experimental results until more experimental data are 
available for comparison. This is due to the fact that the upstream data 
were not available for Bradley [9] and the flow field can be significantly 
affected by the upstream turbulence level, the magnitude of surface 
roughness and its distribution, and the geometry of the obstacle. This 
is demonstrated by the present study. 
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Figure 6.21 indicates that the distribution of turbulence 
kinetic energy is significantly affected by the aspect ratio of the 
obstacle, as well as by the magnitude of the surface roughness scale. 

As mentioned earlier, the length scale, L, on top of the obstacle is 
under-predicted by a prescribed mixing length hypothesis. An under- 
prediction of L causes an under-prediction of the production rate of 
turbulence kinetic energy, K, and an over-prediction of the dissipation 
rate of K, Equation 3.9. Therefore, the results of Taylor [39] show 
relatively lower turbulence kinetic energy levels than the present 
computed results. However, his results agree well with the data of 
Bradley [9]. The two-dimensional model of Taylor [39] developed for 
flow over gentle terrain results in good agreement with the K/Kq profile 
for a set of three-dimensional, full-scale field data for a high hill. 
Figure 6.20, and in disagreement with the V/V^ profile, Figure 6.19. 
Therefore, more data are needed to verify the accuracy of Taylor's 
model when it is applied to the cases of flow over high hills. 

The present results show maximum values of turbulence kinetic 
energy occur at z = 0.3. Taylor's method [39], however, shows a 
constant turbulence layer due to his assumption of a constant shear 
layer in the wall region. His results [39] shown in Figure 6.20 were 
obtained for the case Zg = 0.0005. Since Kg = (V*/C^) is used in the 
present study, Kg = 0.016 if = 0.416 and V* = 0.053, which corresponds 
to the case Zg = 0.0005. The present computed results show that the 
value of K can be as large as ten times its undisturbed value. Kg, 
if Zg = 0.0005 and D/H = 2.0, Figure 6.21. Therefore, the maximum 
value of K is 0.16 if Kg = 0.016. The square root of K = 0.16 results 
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in the value of the ratio »^/V^ = 0.4. Some evidence [5,7] exists 
that the value of Oy^/V^ is approximately 0.4 on top of surface 
obstacles. As discussed earlier, the value of is 1.03 for flow 

over uniform, level terrain. Equation 6.3, and is 0.93 in the layer 
directly on top of Bradley's hill [9]. The value of = 0.4 

obtained in the present computation is reasonable if the value of 
Oy^/V^ = 0.4, which results in the ratio Oy^/ZK = 1.0. 

6.5 Phenomenon of Flow Separation 

As flows pass over a gradually sloping surface, fluid elements 
close to the surface gradually lose momentum due to surface friction. 

Flow separation occurs when the momentum of the fluid particle cannot 
overcome the adverse pressure gradient caused by the obstruction to the 
flow, i.e., the fluid elements flow in the direction of decreasing pres- 
sure. As contrasted to the gradually sloping surface, an abrupt or 
bluff-shaped obstacle creates a small separation on the front face due to 
the same phenomenon described above, i.e., the interaction of adverse 
pressure gradients and the viscous forces; however, a much larger 
downstream separation emanates from the sharp leading edge due to the 
inability of the flow to negotiate the change of surface configuration 
[1]. For flow over non-bluff bodies, e.g. , flow over semiel 1 ipses , 
the mechanism of flow separation and reattachment is the same as that 
of flow over gradually sloping surface [1]. The flow parameters 
affecting the phenomenon of flow separation are therefore the scale 
of the viscous forces, the surface friction, and the geometry of the 
obstruction. In addition to these parameters, the magnitude of the 


upwind turbulence level significantly influences the size of the 
separated flow regions (see Section 6.4). In the present study, 
the above mentioned parameters are scaled by the Reynolds number, 

Re, the surface roughness, Zq, the aspect ratio, D/H, and the 
dimensionless friction velocity, V*, respectively. The effect of 
these parameters on flow separation is discussed in this section. 

The obstacles considered are a rectangular block and different aspect 
ratio semiel 1 iptical cylinders. 

The characteristics of wind around a rectangular block have 
been investigated [11] by using the same numerical model as that used 
in this study. For all the cases studies in Reference [11], the 
presence of a block in the wind field always induces an upstream recir- 
culating flow region and a downstream separation bubble. In general, 
the upstream bubble is small and extends approximately 1 H upstream. 

The size of the downstream separation bubble, however, is significantly 
affected by the upstream turbulence level. In the present model, the 
undisturbed upstream turbulence kinetic energy, Kq, is proportional to 
(V*) . The increase of V* causes a smaller downstream bubble. As 
mentioned in Section 6.3, an increase in Zq of the obstacle tends 
to cause greater displacement of the flow approaching and passing 
over it. Therefore, a larger value of Zq induces a larger upstream 
separation bubble for flow over a semiellipse [19]. However, 
the computation of flow over a forward-facing step [32] shows a minimum 
upstream separation bubble at Zq = 0.05. In view of the fact that the 
quantities obtained in numerical computation show an approximate rather 
than an exact location of the flow separation and reattachment points. 


108 



the effect of Zq on the upstream flow separation should be resolved by 
experimental work. 

Summarily, the flow field upstream of the block is influenced 
by the surface friction, the pressure gradients, the viscous force, and 
the turbulence level of upstream wind. The downstream flow field, how- 
ever, is influenced by the upstream turbulence level and the shape of 
the obstacle. 

The cases of flow over semi elliptical configurations studied in 
this effort did not produce flow separation. The molecular viscosity, 

V, is assumed to be much smaller than the eddy viscosity, therefore, 
the influence of the molecular Reynolds number is not felt. For 

5 

D/H = 1.02, Zq = 0.001, the present computed results with Re = 10 have 

been compared with the case of Re = and only insignificant effects 

2 

on the flow field occur. Assigning conditions Re = 10 , however, 
results in downstream flow separation. Figure 6.22. The flow reattach- 
ment distance, x^, where x^, is measured from the lower corner of the 
ellipse (Figure 6.22) is approximately 2 H, which is much smaller than 
that for laminar flow over a circular cylinder (x^ = 10 H has been 
reported [42]). In view of the fact that the turbulence model described 
in Section 3.0 is developed under the assumption of high Reynolds num- 
bers, reliability of the results obtained by applying the model to low 
Re flows is uncertain. However, assuming an eddy viscosity of zero and 
solving the system of equations described herein, a reattachment distance, 
x^, which extends to 7 H for D/H = 1.0001, is found. This value may 
still be an under-prediction. A coarse grid size used in the computa- 
tion procedure can result in under-prediction of the downstream separa- 
tion bubble of a circular cylinder [42]. 
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Figure 6.23 shows the w contours corresponding to the case 
shown in Figure 6.22. The distribution of the w is very similar to 
the results of a laminar flow calculation for flow over a circular 
cylinder [42]. This implies that the turbulence model developed under 
the assumption of high Reynolds number can be applied to low Reynolds 
number flows by assuming a very small value for turbulence kinetic 
energy, which leads to a small value of It is noteworthy, however, 
that K should not be zero; otherwise a singularity occurs at the term 
of Equation 3.14. 

For laminar flow, increasing Re increases the size of the flow 
separation bubble [42]. An increase in Re can be generated by increas- 
ing the upstream velocity and therefore increasing the momentum of the 
approaching flow. This increased inertia force convects the separation 
bubble downstream. In the case of turbulent flow, an increase in V|^ 
causes a decrease in V*/V^, if V* is kept constant. The downstream 
separation bubble is therefore increased. Figure 6.17. The mechanism 
governing the size of the flow separation bubble, however, is different 
in the laminar flow case and the turbulent flow case. In low Re flows, 
the size of the reversed flow region is controlled by viscous diffusion 
of momentum rather than by turbulent diffusion. The shear layer shed 
behind the body spreads more rapidly and reattaches sooner in turbulent 
flow. 
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Figure 6.23 Vorticity, w, patterns over a semiellipse for D/H = 1.02, 

z =0.001, and Re = 100 
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7.0 CONCLUSIONS 


The values of the constant coefficients developed in the present 
study for use in the K-L turbulence model are applied to the solution of 
the problems of atmospheric flow over two-dimensional obstacles. Good 
agreement between the numerically computed results and the experimental 
results from both atmospheric boundary layer wind tunnels and full-scale 
field tests are obtained. 

The results of studying the influence of the flow parameters 
on the flow field for the range of parameters investigated in this 
study lead to the following conclusions: 

1. The character of the flow field is significantly affected 
by the surface roughness length scale, Zg/H; the upwind 
turbulence level, V*/V^; and the geometries of the obstacles. 

2. The drag force exerted on the flow field by the solid sur- 
face is larger for increased Zg/H and thus more effectively 
retards the flow near the surface, causing the streamlines 
to be further displaced from the obstacle. 

3. The turbulence intensity of the fluid motion is increased in 
the wall region for increased Zg/H. In turn, the momentum 
of the fluid elements in the wall region is reduced due to 
the energy being extracted from the mean flow field and 
converted to turbulence energy. 

4. The thickness of the layer of atmosphere next to the ground 
where the flow is significantly affected by the local surface 
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roughness depends on the distribution of the roughness and 
its magnitude. This thickness is generally less than one 
characteristic body height. 

5. Flow separation generally commences from the upstream corner 
of the bluff body. For flow over a convex surface, the flow 
separates at a position downstream of the crest of the 
obstacle if Re is low. For high Re, the flow is not likely 
to separate, particularly if the upstream turbulence is high. 

6. An increase in upstream turbulence level reduces the size 
of the downstream flow separation bubble. This effect is 
attributed to enhanced turbulent mixing which transports 
momentum much more quickly from the upper layer to the 
lower layer. 

The results of the present study indicate that the turbulent 
flow field upstream of an obstacle is mainly affected by the surface 
friction and the pressure gradient. Downstream of the obstacle, however, 
the flow is mainly affected by the upstream turbulence level. Addition- 
ally, the mean velocity field and the turbulence field are significantly 
affected by the geometry of the obstacle. For low Reynolds number 
flows, the size of the downstream reversed flow region is controlled 
by viscous diffusion of momentum rather than by turbulent diffusion. 
Therefore, the shear layer shed behind the obstacle spreads more 
rapidly in turbulent flow, resulting in a smaller downstream separation 
bubble. 
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APPENDIX 


EXPRESSION OF S^^ IN ORTHOGONAL CURVILINEAR 
COORDINATE SYSTEMS 


Consider an orthogonal curvilinear system (u^u^) rotated an 
angle e with respect to the Cartesian system (x,z) as shown in 
Figure A. 1(a), the unit vectors e and e in (x,z) system can be 

X "Z. 

expressed as: 

e^ = cos 0 e^ - sin 0 e ^ 

e^ = sin 0 e-| + co: 0 e^ (A.l) 

where e^ and are the unit vectors in the (u-ijU^) system. 

The distances between two nearby points in the (u-|,U 2 ) system, 
i.e., dJi and dn in Figure A. 1(b), can be written as 


da = h-|du-j 

dn = h^du^ (A. 2) 


where h^ and h^ are metric coefficients in (u^jU^) system. The deriva- 
tive of a quantity (p with respect to u^ can be expressed as: 


IsL 

1 = fli' 

9x 9£ 

1 + N 

9Z 

3U,J 


dZ 9Ut 

1 ^ 


95, 9U^^ 


= h 


1 


cos 0 


3X 


+ sin 0 


9Z 


(A. 3) 


Similarly 
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z 



(a) Two-dimensional orthogonal curvilinear coordinate 
systems 



Figure A.l Schematic diagram of the relationship between general 
orthogonal curvilinear coordinate system and Cartesian 
coordinate system 


120 



(A.4) 




f-sin 6 


+ cos 0 
8x 


9 (j) 


Multiplying Equations A. 3 and A.4 with cos 0 and h-j 
tively, and then subtracting, one obtains 

9(() _ COS 0 9()) sin 0 9(f) 

9x ~ h-| 9u-| ~ Su^ 

In turn, multiply Equations A. 3 and A.4 with h^ sin 0 
respecti vely, and adding gives 


M = sin 0 9(f) cos 0 9(j) 
9z ~ h-| 9u-| 9U2 


Substituting Equations A.l, A. 5 and A. 6 into Equation 
sion for V and V in the {u-|,U 2 ) system becomes: 


V(|) 


h-| 9u-| h^ 9U2 


V(J) 


-4- 

Al- I ^2 9(j) 
^2 ^^2 ^1 ^^1 


^ • A 

The velocity vector V appearing in Equation 4 


as: 


V = V^e^ + 72^2 


Expanding Equation 4.30 with the aid of Equations A.l 


sin 0, respec- 

(A.5) 

and h^ cos 0, 

(A. 6) 

4.31, the expres- 
(A.7) 

.30 can be written 
(A. 8) 

, A. 7 and A. 8 
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9y 

9V^ 9y^ 


X 

X J X 

X 

9 U 2 

9 U-j 9 u-j 

9 U 2 

9y^ 

aV^ 9y^ 

9V 

z 

^ + ^ 

z 

3U2 

9U-| 9U-j 

9 U 2 


where 


9u . 9u 

COS e e _ sin e e 

h-| 9u-| h2 


• 9y 

sin 9 e cos Q e 

h-| 9U-J 9 u2 

V = Vt cos 0 - Vo sin 0 

X 1 2 

V = V, sin 0 + Vo cos 0 

z 1 2 

Substituting Equation A. 10 into A. 9 and arranging gives Equation 
the text: 

S 9V, 9y, 9V, 9y, 9Vo SUo 9Vo SUo 

hh — = 1—1 + 1 + 

1 2 2 9u-| 9 U 2 9U2 9u-| 9u-| 9u 2 9^2 3u-j 



r 

90 

3V2 

90 

9 V 2 ' 

+ ^2 

90 

9V^ 

90 

3VJ 

1 

9U^ 

3U.2 

9u2 

9Ui 

' j 

902 

< 

9u-j 

9u-| 

902 

^ J 


90 

3y2 

90 

302 






1 

9 U 2 

s. 

9U-| 

9u-| 

9U2 

j 







90 

9y-| 

90 

9y-j 






2 

9U-| 

8U2 

902 

9u-j 







where 


(A.9) 


(A. 10) 
.32 in 


(A. 11) 
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= 

^2 fi2 3U2 


(A. 12) 
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